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A  Project  Abstract 


Frequency  hopping  (FH)  is  the  prevailing  spread  spectrum  method  in  military  communications, 
largely  due  to  its  low  probability  of  detection  and  interception.  In  this  project  we  developed  a  novel 
signal  processing  scheme  for  code-blind  reception  of  multiple  frequency  hopped  transmissions 
over  multipath  channels.  This  technique  is  based  on  the  principle  of  dynamic  programming  and/or 
expectation-maximization,  coupled  with  multidimensional  harmonic  and  low-rank  analysis.  It  is 
able  to  jointly  estimate  hop  timing,  hop  frequency,  and  direction-of-arrival  (DOA)  of  multiple  FH 
signals  in  the  presence  of  frequency  collisions,  without  the  knowledge  of  signal  hop  patterns.  The 
method  can  also  be  used  to  obtain  locate  information  and  operation  characteristics  of  active  FH 
jammers.  The  associated  identifiability  of  2-D  and  multidimensional  frequency  estimation  is  also 
investigated. 


B  Technical  Report 

B.l  Problem  Statement  and  Summary  of  Major  Results 

This  project  studied  the  problem  of  detection,  localization,  and  tracking  of  multiple  frequency 
hopped  signals  in  multipath  channels,  without  the  knowledge  of  their  hop  codes  and  hop  timing. 
Interception  and  localization  of  multiple  FH  signals  are  challenging  problems  that  feed  into  mul¬ 
tiple  facets  of  military  communications,  from  interception  of  noncooperative  communications  to 
jammer  mitigation.  On  the  technical  side,  the  problems  are  challenging  not  only  because  hopping 
patterns  and  DOAs  are  unknown,  but  also  other  parameters  such  as  frequency  bin- width,  hop-rate, 
and  timing  are  at  least  partially  unknown  in  a  realistic  scenario.  Carrier  hopping  means  that  one 
has  to  deal  with  switching  exponentials,  rather  than  pure  exponentials;  and  it  also  induces  hopping 
in  the  receive  antenna  array  spatial  steering  vectors,  due  to  wavelength-dependent  phase  shifting 
from  one  array  element  to  another. 

Many  signal  processing  techniques  have  been  developed  for  blind  interference  suppression  in 
frequency  hopping  using  an  antenna  array,  e.g.,  [26,29-31].  These  approaches  aim  for  interference 
suppression  rather  than  joint  multiuser  detection  and  hopping  pattern  identification,  they  all  require 
at  least  knowledge  of  the  hopping  pattern  of  a  given  signal  of  interest,  and  their  interference  nulling 
capability  is  bounded  by  the  degrees  of  freedom  in  the  adaptive  array.  On  the  other  hand,  several 
papers  have  been  published  on  the  subject  of  (joint)  multiuser  detection  for  frequency  hopping 
systems,  e.g.,  [4, 18].  These  assume,  among  other  things,  that  the  hopping  patterns  of  all  signals 
are  known  to  the  receiver,  hence  clearly  not  applicable  in  a  noncooperative  scenario. 

Without  assuming  the  knowledge  of  hop  patterns,  several  methods  have  been  proposed  for 
blind/semi-blind  hop  timing  and  frequency  estimation.  For  example,  assuming  known  hop  rate, 
channelized  receivers  have  been  proposed  for  semi-blind  hop  timing  estimation  (knowledge  of 
frequency  channelization  is  required)  for  the  single  user  case  [2,  25],  as  well  as  the  multiuser 
case  [1],  However,  the  performance  of  those  receivers  degrades  rapidly  if  the  channelization  is  im¬ 
perfect.  [32]  considers  adaptive  source  localization  and  blind  beamforming  for  FH  signals  without 
the  knowledge  of  hopping  patterns.  The  algorithm  requires  rough  synchronization  with  the  desired 
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signal’s  hop  interval  and  is  restricted  to  the  use  of  very  special  array  geometries,  and  the  number  of 
users  that  can  be  resolved  is  bounded  by  the  number  of  sensor  array  elements.  In  our  preliminary 
work  [17],  a  first  step  is  used  to  separate  multiple  FH  signals  and  then  a  single  user  is  tracked  based 
on  1-D  harmonic  retrieval. 

The  major  results  of  this  project  are  summarized  as  follows. 

1.  We  developed  the  DP-2DHR  algorithm,  based  on  the  principle  of  dynamic  programming 
(DP)  coupled  with  2-D  harmonic  retrieval  (2-D  HR),  to  jointly  estimate  the  hop  timing, 
hop  frequency,  DOA,  of  multiple  FH  signals  in  multipath  channels  with  possible  frequency 
collisions.  The  algorithm  does  not  require  the  knowledge  of  hop  patterns  or  channelization. 

2.  Using  a  multiple  invariance  sensor  array,  we  developed  the  DP-TALS  algorithm  with  low- 
rank  decomposition  of  three-dimensional  arrays,  for  joint  timing,  frequency,  and  2-D  DOA 
estimation  of  multiple  FH  signals  without  hop  pattern  knowledge. 

3.  Relying  the  principle  of  expectation-maximization  (EM),  we  developed  a  low-complexity 
EM  algorithm  for  joint  hop  timing  and  frequency  estimation  by  exploiting  the  inherent  data 
structure  of  multiple  frequency  hopped  transmission. 

4.  We  advanced  the  theory  of  identifiability  of  2-D  and  multidimensional  frequency  estimation, 
and  proved  the  most  relaxed  identifiability  bound  of  N-D  frequency  estimation  to  date.  An 
eigenvector-based  algebraic  algorithm  for  2-D  frequency  estimation  was  also  developed, 
which  achieves  the  identifiability  bound,  and  offers  asymptotically  optimal  performance. 

B.2  The  DP-2DHR  Algorithm  for  Timing  and  Frequency  Estimation 

B.2.1  Data  Modeling 

Suppose  an  M-element  uniform  linear  array  (ULA)  receives  FH  transmission  from  d  sources.  Each 
far  field  FH  signal  is  from  a  nominal  DOA  with  negligible  angel  spread.  The  array  steering  vector 
in  response  to  a  signal  from  direction  a  is  a (6)  =  [1,  6, ... ,  9M~1]T,  where  9  =  ei27rAsm(«).  The 
received  signal  is  sampled  at  a  sampling  rate  of  1/T  ( T  is  normalized  to  1),  and  the  M  x  1  signal 
vector  collected  at  the  ULA  output  at  sampling  time  n  can  be  expressed  as 

d 

x{n )  =  ^2  a(^p))/3^p)sr(n)  +  w(n ),  (1) 

r— 1 


where  sr(n )  =  eJuj(,P>r\  and  is  the  frequency  of  the  transmitted  signal  from  the  r-th  user  during 
its  p-th  hop.  Note  that  the  baseline  separation  A  (measured  in  wavelength  units)  is  frequency 
dependent,  hence  so  are  the  steering  vectors.  For  notational  clarity,  sometimes  we  do  not  explicitly 
denote  this  dependence  as  long  as  it  is  clear  from  the  context.  The  transmitted  signals  can  be 
fast  or  slow  frequency  hopping,  with  FSK  or  linear  modulation.  3r'}  is  the  complex  path  loss 
for  the  r-th  user  during  its  p-th  hop  that  collects  the  (frequency-dependent)  channel  attenuation; 
the  signal’s  initial  phase  (jhl>}  is  also  absorbed  into  diP! .  Here  the  carrier  shifts  due  to  hopping  or 
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Figure  1:  An  example  of  two  frequency-hopped  signals. 

symbol  modulation  are  treated  as  conceptually  equivalent,  albeit  of  different  magnitude.  w(n)  is 
complex  white  Gaussian  noise  with  variance  a2.  Suppose  N  samples  (snapshots)  are  collected  at 
the  array  output,  then  the  received  data  matrix  can  be  written  as  X  =  [a?(0)  ■  ■  ■  x(N  —  1)]. 

The  objective  here  is  to  estimate  hop  timing  (i.e.,  hop  instants)  and  hop  frequency  sequences 
of  all  transmitted  signals  from  X,  in  the  presence  of  possible  collisions  in  some  time  segments, 
without  the  knowledge  of  signals’  hop  patterns  or  hop  rates.  In  model  (1)  we  assume  a  single-path 
transmitter-receiver  propagation  for  each  signal.  Later  we  generalize  it  to  incorporate  multipath 
propagation  with  small  delay  spread.  We  also  assume  that  the  number  of  signals  (and,  in  a  multi- 
path  scenario,  the  total  number  of  paths  for  all  signals)  has  already  been  estimated  by  an  appropriate 
source  enumeration  method,  such  as  rank  criteria  (e.g.,  SVD)  or  information  theoretic  criteria  (e.g., 
AIC,  MDL). 

For  simplicity  of  exposition,  let  us  focus  on  an  FH  system  where  there  are  two  active  FH 
signals.  As  shown  in  Figure  1,  si(t)  and  s2(^)  may  have  different  hop  rates  and  hop  timing.  Let 
nt,  i  =  0, . . . ,  K  —  1,  be  the  system-wide  hop  instants  (n0  =  0  and  uk  =  N  by  convention).  We 
assume  that  within  time  period  of  interest,  the  total  number  of  hops  of  all  signals  is  bounded  above 
by  K  —  1  (such  a  bound  could  be  deduced  from  the  spectrogram  of  the  data,  and  need  not  be  tight). 

Between  any  two  system-wide  consecutive  hop  instants,  e.g.,  n7  and  ni+ 1,  there  are  only  two 
temporal  frequencies  involved.  During  such  a  time  segment,  the  received  data  may  be  written  as 

X i  =  [; x(rii )  •  •  •  x(ni+ 1  -  1)]  =  AiBiSj  +  Wh  (2) 

where  A,  =  \a(f)ip})  a{d^)\,  B,  =  diag (H^),  and  the  subscript is  a  time  index  indicating 
that  the  time  segment  is  delimited  between  nt  and  nl+  \  —  1,  i.e.,  the  z-th  system-wide  dwell.  In  (2), 
the  signal  matrix  Si  is  defined  as 

eju>[p\ni+l)  .  .  .  ejuj[p\ni+ i  —  l)  T 
ejoj^\rii+l)  .  .  .  eju(2q\ni+  i  —  l) 

and  W i  is  the  corresponding  noise  matrix.  Here  we  assume  user  1  and  user  2  are  in  their  p-th 
and  g-th  hops  respectively  during  this  time  segment,  and  a(()\pl)  and  a (O^)  are  the  antenna  steer¬ 
ing  vectors  corresponding  to  and  Since  both  Aj  and  S,  are  Vandermonde  matrices, 

the  estimation  of  DOAs  and  frequencies  from  X,  in  (2)  is  in  fact  a  2-D  constant  modulus  har¬ 
monic  retrieval  problem,  and  there  are  two  frequency  components  along  each  of  the  spatial  and 
temporal  dimensions.  If  d  users  are  active  in  the  system,  a  similar  2-D  harmonic  mixture  model 


S7  = 


■  (p) 

rii 

■  ( q ) 

gJw2  ni 
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can  be  obtained  except  that  the  number  of  frequency  components  in  such  a  time  segment  along 
each  dimension  is  d.  Recently,  improved  identifiability  results  and  algorithms  regarding  2-D  HR 
have  been  developed  [6, 8, 11, 14, 15].  We  use  the  MDF  algorithm  in  [15]  for  the  purpose  of  2-D 
frequency  estimation  here. 

B.2.2  Joint  Timing  and  Frequency  Estimation 

The  key  idea  behind  our  proposed  method  of  joint  hop  timing  and  frequency  estimation  is  that 
between  any  two  hypothesized  system-wide  hops,  the  data  follow  a  2-D  harmonic  model.  Hence 
for  a  hypothesized  set  of  hops  (that  is,  including  all  hops  of  all  users  in  the  system),  2-D  harmonic 
retrieval  methods  can  be  used  to  estimate  model  parameters,  and  subsequently  calculate  model 
fit.  If  one  operates  under  an  upper  bound  on  the  total  (system- wide)  number  of  hops,  then  system 
stage  can  be  defined  as  the  number  of  allowable  hops,  and  state  can  be  defined  as  the  hop  instant, 
hence  dynamic  programming  can  be  used  to  find  the  optimal  hop  sequence  and  associated  model 
parameters  per  dwell  [16].  With  d  users  and  a  budget  of  K  —  1  hops,  define 


n  =  [nt, . . .  ,nA-_i], 


OL  = 

[a?",- 

(K- 1) 
ry )  J 

.  .  ,  9  • 

(0) 

..,ayd\. 

(K-  lf- 
ry  \ 

•  '  } 

P  = 

[/sf1,. 

6(0) 

■  ■  iPd  i- 

<jl!  = 

M0),. 

,  ,(0) 

■  ■  iud  v ' 

as  the  vectors  of  hop  timing,  DOAs,  complex  frequency-dependent  attenuations,  and  hop  frequen¬ 
cies.  Joint  maximum  likelihood  estimation  of  n,  ck,  / 3 ,  and  u>  from  X  amounts  to  minimizing 

K—l 

J(h,  a,  /3,  u?)  =  ||Xj  -  Xj|||  (3) 

i= 0 

over  n,ct,(3,i2>,  where  X,  is  the  reconstructed  2-D  harmonic  mixture  based  on  ML  parameter 
estimates  (DOAs,  complex  amplitudes,  and  carrier  frequencies),  obtained  in  each  time  segment 
defined  by  hypothesized  hl  and  nl+  \ ,  assuming  a  2-D  harmonic  mixture  model  for  the  received 
data  during  this  segment.  Since  K  will  typically  be  higher  than  the  true  number  of  hops  in  the 
available  samples,  we  include  a  “parking  stage”  in  the  DP  program  to  account  for  the  possibility 
of  unused  hops.  In  the  presence  of  noise,  however,  DP  will  typically  use  any  extra  hops  available 
to  track  minor  noise-induced  variations.  Such  variations  can  be  relatively  easily  detected  after  DP, 
for  frequencies  before  and  after  such  hops  will  be  approximately  equal. 

From  the  MDF  estimates,  we  form 

d 

x(n )  =  J2a(8r)Prej*rn, 

r= 1 

for  rii  <  n  <  n,+ 1 ;  here,  6r  =  ej27rAsm(“r),  and  the  matrix  X,  is  constructed  from  x(n)  in  the  same 
form  as  Xj  in  (2).  Define  nl+i  —  1],  for  0  <  i  <  K  —  1,  as  the  cost  function  for  the  time 
segment  n,  <  n  <  ni+ 1 

Aj[rii,  ni+i  —  1]  =  \\Xi  —  Xi\\F.  (4) 
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Table  1 :  The  DP-2DHR  Algorithm 


1.  Initialization 

Let  k  =  1,  compute  Tk(L)  for  L  =  2, . . . ,  N  —  3 K  +  2  using  Eqns.  (4) 
and  (5). 

2.  Recursion 

For  2  <  k  <  K  —  1,  compute  Tk(L)  with  L  =  3k  —  1, N  —  3 K  + 
3k  —  1,  using  Eqn.  (6)  with  3k  —  3  <  nk~ 1  <  L  —  2\  For  A;  =  A\  compute 
Tfc(L)  with  L  —  N  —  1. 

For  each  L,  denote  the  value  of  that  minimizes  T k ( L )  as  /</,•- 1  (L), 

and  denote  the  corresponding  «*._!,  as  ak_i(L), 

and  u>k_i(L),  respectively. 

3.  Backtracking 

The  maximum  likelihood  estimates  of  hop  instants  are  obtained  by  using 
the  backward  recursion,  i.e.,  hi  =  rii{hi+ 1  —  1),  for  i  —  K  —  2,  K  — 
3, . . . ,  1,  initialized  by  fiK- i  =  nx- i(N  —  1).  Similarly,  the  corre¬ 
sponding  DOA,  amplitude,  and  frequency  estimates  of  each  segment  can 
be  obtained  by  their  respective  backward  recursions. 


Furthermore,  to  solve  the  minimization  problem  in  (3)  by  DP,  we  define 

fc-i 

Tk(L)=  min  V'  AJrij,  ni+1  -  1],  (5) 

nl  —  1  L ^ 

n0=0,nk=L+l  i=0 

where  0  <  rii  <  •  •  •  <  nk-\  <  L.  Eqn.  (5)  can  be  viewed  as  the  minimization  problem  of  finding 
the  best  fit  for  a  subset  of  the  data  of  size  M  x  (L  + 1)  when  a  total  number  of  k  —  1  hops  is  allowed. 
Hence  Tk(N  —  1)  is  the  minimum  of  J(h,  a,  u ),  0).  From  (5),  a  recursion  for  the  minimum  can 
be  developed  as 

Tfc(L)  =  min  -  1)  +  Afc_i[nfc_i,  L] ) .  (6) 

This  means  that  for  a  data  matrix  of  size  M  x(L  +  1),  the  minimum  error  for  k  segments  (i.e.,  k  —  1 
hop  instants)  is  the  minimum  error  for  the  first  k  —  1  segments  that  end  at  n  =  nk_ \  —  1,  and  the 
error  contributed  by  the  last  segment  from  n  =  nk_\  to  n  =  L.  The  solution  of  the  minimization 
of  (3)  is  for  k  =  K  and  L  —  N  —  1,  which  yields  the  joint  estimates  of  hop  timing,  DOAs,  and 
frequencies  of  all  signals. 

Assuming  that  the  minimum  length  of  a  segment  is  three  samples,  the  procedure  to  compute 
the  solution  by  the  DP  and  2-D  Harmonic  Retrieval  (DP-2DHR)  algorithm  is  summarized  in  Table 
1.  Note  that  frequencies  and  complex  amplitudes  of  different  segments  pertaining  to  a  particular 
signal  can  be  associated  via  their  corresponding  DOA  parameters,  since  for  a  single  segment, 
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frequency  and  DOA  parameters  pertaining  to  one  signal  are  paired  up  automatically  by  the  MDF 
algorithm.  Depending  on  different  transmission  schemes,  the  application  of  the  DP-2DHR  method 
may  slightly  vary  as  described  in  the  following  cases. 

1.  Slow  frequency  hopping  (SFH)  with  M-ary  FSK  modulation:  Frequency  changes  due  to 
baseband  modulation  are  usually  much  smaller  than  those  due  to  carrier  hopping.  Hence 
symbol  rate  and  hop  rate  can  be  obtained  from  the  result  of  DP,  and  consequently  symbol 
recovery  is  possible. 

2.  SFH  with  M-ary  PSK  or  M-ary  QAM  modulation:  During  one  hop  dwell,  frequency  is 
constant,  but  the  complex  amplitudes  are  different  from  symbol  to  symbol  due  to  modulation 
(recall  that  for  one  hop  dwell,  the  effect  of  channel  on  the  complex  amplitudes  is  constant). 
Hop  timing  can  be  detected  from  frequency  change.  Hence  symbol  rate  and  hop  rate  are 
distinguishable  from  the  result  of  DP. 

3.  SFH  with  GMSK  modulation:  A  GMSK  signal  is  not  a  pure  exponential  in  one  symbol 
period.  However,  narrowband  GMSK  can  be  well- approximated  by  a  pure  exponential  for 
our  purpose. 

4.  Fast  frequency  hopping  (FFH):  The  DP-2DHR  method  is  applicable  for  hop  timing  and 
hop  frequency  sequence  estimation.  However,  additional  information  is  needed  for  symbol 
detection,  e.g.,  symbol  period  and  symbol  synchronization  are  required  since  the  DP-2DHR 
can  only  provide  chip  synchronization  in  this  case. 

Next  we  present  the  numerical  simulation  results  to  demonstrate  the  proposed  DP-2DHR  algo¬ 
rithm  for  joint  hop  timing  and  frequency  estimation  in  the  presence  of  frequency  collisions.  Two 
FH  signals  with  DOAs  [12°,  17°]  are  simulated,  each  hopping  with  different  hop  timing.  The  re¬ 
ceiver  array  has  M  =  6  antennas,  with  baseline  separation  of  A/2  at  fc  =  1  GHz.  With  M  =  6, 
the  array  has  a  3dB  beamwidth  of  about  28°  so  that  the  two  sources,  separated  by  5°,  are  not  di¬ 
rectly  resolvable.  A  hopping  frequency  band  of  bandwidth  8  MHz  is  occupied  by  32  frequency 
channels  with  0.25  MHz  channel  spacing.  The  received  signal  is  well-modeled  as  narrow-band. 
For  simplicity  of  illustration,  hop  rate  is  set  the  same  as  symbol  rate  (125  Kbps).  At  the  receiver, 
the  complex  antenna  outputs  are  sampled  at  a  rate  of  8  MHz  after  down-conversion,  and  N  =  48 
complex  samples  are  collected  at  each  antenna,  resulting  in  a  6  /is  long  analysis  window,  hence 
each  signal  hops  at  most  once  within  this  window.  SNR  is  defined  as 

SNR:=ioiog«(!^l)’  m 

where  the  noise  variance  a2  =  N0B,  and  B  is  the  processing  signal  bandwidth. 

An  example  for  the  FH-FSK  case  is  shown  in  Table  2  and  Figure  2.  In  this  example,  two  BFSK 
signals  begin  in  different  bins;  then  signal  1  hops  to  the  same  bin  as  signal  2;  then  signal  2  hops 
out  of  his  original  bin  and  into  a  new  bin.  This  gives  three  segments:  the  first  and  the  third  without 
collisions,  and  the  second  with  collisions.  Table  2  gives  the  DOA  estimation  results  for  the  three 
segments  and  Figure  2  gives  the  corresponding  results  of  hop  timing  and  frequency  estimation 


8 


Table  2:  DP-2DHR:  DOA  estimation  of  two  FH-FSK  signals. 


True 

DOA 

Estimated  DOA 

1st  Seg. 

2nd  Seg. 

3rd  Seg. 

Signal  1 

12° 

12.24° 

13.23° 

12.69° 

Signal  2 

17° 

17.35° 

16.10° 

17.16° 

Signal  1 

-  Estimated 

-  -  Original 

. - „„ 

0i - i - ^ ^ - i - i - 1 

0  1  2  3  !  ,  ,  v4  5  6 

:  t  (/is) 
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Figure  2:  DP-2DHR:  multiuser  tracking  of  FH-FSK  signals  (SNR=10dB). 

for  the  two  signals.  SNR  is  lOdB.  We  assume  that  any  mobility-induced  changes  in  DOA  are 
negligible  within  the  analysis  window.  Thus,  varying  hop  frequencies  are  associated  with  different 
signals  via  their  corresponding  window  invariant  DOA  parameters.  The  results  show  that  DOA, 
hop  timing  and  frequency  estimates  are  close  to  the  respective  true  values  even  in  the  presence  of 
collisions.  They  also  demonstrate  that  good  estimates  can  be  obtained,  based  on  measurements  of 
duration  less  than  one  symbol  period;  this  implies  that  the  algorithm  is  capable  of  blind  multiuser 
tracking  at  moderate  to  heavy  loads. 

Figure  3  depicts  the  Root  Mean  Square  Error  (RMSE)  of  DP-2DHR  hop  timing  and  frequency 
estimates  in  the  presence  of  collisions.  The  RMSE  is  obtained  via  Monte  Carlo  simulation.  For 
each  realization,  each  of  the  two  FH-FSK  signals  hops  once  within  the  observation  window.  Hop 
timing  is  randomly  generated,  and  frequencies  are  also  randomly  selected  from  the  32  candidate 
bins  with  the  constraint  that  there  is  always  one  collision  in  the  three  hop-free  segments.  The 
RMSE  vs.  SNR  performance  shown  in  Figure  3  indicates  that  the  DP-2DHR  algorithm  performs 
quite  well  even  in  the  low  SNR  regime,  given  the  fact  that  the  signals  are  tracked  in  a  situation 
where  hop  pattern,  rate,  and  timing  are  all  unknown. 

B.2.3  Multipath  Channels  with  Small  Delay  Spread 

The  DP-2DHR  method  developed  in  Section  B.2.2  assumes  single -path  transmitter-receiver  prop¬ 
agation  for  each  frequency  hopped  signal.  When  the  signal  bandwidth  is  greater  than  the  channel 
coherence  bandwidth,  channel  effects  due  to  multipath  propagation  cannot  be  ignored.  Multipath 
reflections  create  fictitious  sources  in  the  spatial  dimension,  as  well  as  unknown  delay  spread  in 
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Figure  3:  DP-2DHR:  RMSE  of  timing  and  frequency  estimation  vs.  SNR. 

the  temporal  dimension.  Few  techniques  can  be  found  in  the  literature  on  blind  parameter  estima¬ 
tion  for  FH  signals  in  multipath  channels,  e.g.,  [7],  which  assumes  a  fixed  but  unknown  header  in 
each  packet  of  the  FH  signal.  In  this  section  we  show  that  the  DP-2DHR  method  can  be  extended 
to  blindly  estimate  hop  timing  and  frequency  of  multiple  FH  signals  in  multipath  channels  with 
unknown  but  small  delay  spread. 

Suppose  the  r-th  signal  arrives  at  the  ULA  from  Lr  distinct  paths  due  to  multipath  propagation, 
then  (1)  can  be  extended  to 


d  Ly* 

x(n)  =  a  i9^)  Prfsr  (n  -  Tri )  +  w(n)  (8) 

r=  1  1=1 

=  EE°  («“)  +  Mn)  (9) 

r= 1  /=1 

for  n  —  0, . . . ,  N  —  1,  where  =  (3^ e^P) Trl .  Here  we  assume  that  the  delay  spread  for  a  given 
signal  is  small  so  that  time  delay  can  be  approximated  by  phase  shift. 

Between  any  two  system-wide  consecutive  hop  instants,  e.g.,  nt  and  ni+ 1,  the  received  data 
can  be  expressed  in  matrix  form  X%  =  AiB.Sf  +  W *,  which  is  a  2-D  harmonic  mixture  model 
that  is  essentially  the  same  as  (2),  with  the  difference  that  the  number  of  frequency  components  in 
such  a  time  segment  along  each  dimension  is  L\  +  L-2  +  •  •  •  +  Lj.  Some  frequency  components 
in  the  time  dimension  are  identical  due  to  multipath  reflection,  but  they  can  be  dealt  with  by  the 
MDF  algorithm.  Hence  DP-2DHR  can  be  applied  as  before  for  joint  hop  timing  and  frequency 
estimation. 

Note  that  frequencies  and  complex  amplitudes  of  different  segments  pertaining  to  a  particular 
path  can  be  associated  via  their  corresponding  DOA  parameters,  since  for  a  single  segment,  fre¬ 
quency,  amplitude,  and  DOA  parameters  pertaining  to  one  path  are  paired  up  automatically  by  the 
MDF  algorithm.  In  addition,  different  paths  pertaining  to  a  particular  emitter  will  result  in  differ¬ 
ent  DOAs  but  identical  hop  frequency  sequence  and  hop  timing  (recall  that  time  delay  is  treated  as 
phase  shift),  hence  paths  can  be  associated  with  emitters  by  hop  sequences,  which  is  a  clustering 
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Table  3:  DP-2DHR:  DOA  estimation  in  the  presence  of  multipath. 


True 

DOA 

Estimated  DOA 

1st  Seg. 

2nd  Seg. 

3rd  Seg. 

Path  1-1 

6° 

4.67° 

7.48° 

6.33° 

Path  1-2 

14° 

12.25° 

14.85° 

13.55° 

Path  2-1 

25° 

25.50° 

24.91° 

24.17° 

(a)  Path  1  of  Signal  1 

-  Estimated 

-  -  Original 
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0  1  2  3  4  5  6 
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Figure  4:  DP-2DHR:  multiuser  tracking  in  the  presence  of  multipath. 

problem  and  can  be  solved,  e.g.,  by  calculating  the  pair-wise  distance  among  all  recovered  hop 
sequences. 

An  example  of  multiuser  tracking  by  the  DP-2DHR  algorithm  is  shown  in  Table  3  and  Figure  4. 
Again,  varying  hop  frequencies  are  associated  with  different  paths  via  their  corresponding  window 
invariant  DOA  parameters.  The  results  show  that  hop  timing  and  frequency  estimates  are  close  to 
their  respective  true  values.  Figure  4  also  indicates  that  paths  1  and  2  pertain  to  the  same  emitter 
since  they  have  essentially  the  same  hop  timing  and  frequency  sequence.  Similar  results  have  been 
obtained  for  GMSK  modulated  signals  [16]. 

Figure  5  plots  the  RMSE  of  hop  timing  and  frequency  estimation  of  the  DP-2DHR  algorithm  in 
the  presence  of  multipath  propagation.  The  results  indicate  that  the  DP-2DHR  algorithm  performs 
well  in  multipath  channels,  given  the  fact  that  the  signals  are  tracked  in  a  situation  where  hop  code, 
rate,  timing,  and  multipath  delay  are  all  unknown. 

In  practical  systems,  the  number  of  signals  is  much  less  than  the  number  of  time  samples, 
and  the  number  of  antenna  elements  usually  ranges  from  3  to  8.  It  can  be  shown  that  a  good 
estimate  of  the  complexity  of  the  DP-2DHR  algorithm  is  0(KN5).  There  are  several  ways  that 
this  complexity  can  be  reduced:  i)  It  is  only  during  the  initial  acquisition  period  that  the  full 
complexity  of  the  blind  algorithm  is  needed.  If  frequencies  hop  at  a  regular  rate,  hop  timing  and 
hop  period  can  be  estimated  by  applying  the  DP-2DHR  algorithm  to  a  relatively  short  portion 
of  a  long  data  record,  while  frequency  estimation  for  the  remaining  data  can  be  accomplished 
by  applying  the  MDF  algorithm  to  pre-decided  hop-free  data  blocks  delimited  by  system-wide 
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Figure  5:  DP-2DHR:  RMSE  of  timing  and  frequency  estimation  in  the  presence  of  multipath. 

adjacent  hop  instants.  This  will  reduce  the  complexity  significantly,  ii)  The  DP-2DHR  algorithm 
may  be  simplified  using  standard  approaches  of  reduced-complexity  Viterbi  decoding,  such  as 
path  pruning  based  on  metric  thresholding,  early  path  merging,  etc.  These  will  of  course  incur 
a  performance  loss,  but  if,  e.g.,  the  truncation  parameter  is  appropriately  chosen,  the  loss  will 
be  small,  iii)  If  one  has  a  reasonably  good  idea  about  the  hop  rates,  the  problem  can  be  much 
simplified.  If  the  hop  code  and  hop  timing  of  a  signal-of-interest  are  known,  then  one  can  de-hop 
and  obtain  a  model  with  much  reduced  noise  and  interference,  since  only  interferers  who  collided 
with  the  particular  user  of  interest  within  the  observation  interval  will  remain  in  the  de-hopped 
signal;  and  the  receiver  can  cut-down  its  bandwidth  to  the  hopping  bin-bandwidth  in  this  case,  iv) 
As  a  further  alternative,  simpler  frequency  estimation  techniques  can  be  used  in  place  of  2-D  MDF. 
Clearly,  there  are  many  trade-offs  one  may  pursue. 

In  the  development  of  the  DP-2DHR  algorithm,  signal  bandwidth  is  assumed  to  be  known.  In  a 
practical  blind  estimation  scenario,  the  receiver  may  also  lack  knowledge  of  the  signal  bandwidth. 
Relative  to  the  other  unknowns  (hop  patterns,  timing  and  rates),  it  is  simpler  for  the  receiver  to  es¬ 
timate  the  compound  signal  bandwidth,  e.g.,  via  energy  detection.  However,  due  to  sampling  rate 
and  noise  power  considerations,  an  intercept  receiver  may  only  observe  part  of  the  spread  band¬ 
width.  In  this  case,  the  performance  of  DP-2DHR  will  be  degraded  due  to  bandwidth  mismatch 
because  signals  may  hop  in  and  out  of  the  observed  band,  making  it  difficult  to  track  across  hops. 
Identifiability  issues  also  become  much  more  complicated  in  the  presence  of  bandwidth  mismatch, 
due  to  model  order  variations. 

We  have  conducted  extensive  simulation  experiments  to  study  the  performance  of  the  DP- 
2DHR  algorithm.  More  results  may  be  found  in  [16].  In  addition,  based  on  a  software  testbed 
developed  for  [17],  we  have  designed  a  new  software  demo  for  blind  multiuser  tracking  in  FH 
systems  using  Matlab.  A  screen  shot  of  the  demo  is  shown  in  Figure  6. 

B.3  The  DP-TALS  Algorithm  with  Multiple  Invariance  Sensor  Arrays 

The  principle  of  DP-2DHR  can  be  extended  to  jointly  estimate  2-D  DOA  (azimuth  and  elevation 
angles),  hop  timing,  and  frequency  of  multiple  FH  transmissions  using  an  antenna  array  possess- 
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Figure  6:  A  software  demo  for  blind  multiuser  tracking  in  frequency  hopping  systems. 

ing  multiple  invariances  (MI).  An  MI  sensor  array  is  composed  of  multiple  identical  subarrays 
displaced  in  the  same  or  different  directions.  A  multiple  invariant  rectangular  array  is  shown  in 
Figure  7.  Several  methods  have  been  proposed  for  direction  finding  and/or  other  parameter  esti¬ 
mation  using  MI  sensor  arrays,  e.g.,  [24,27,28,34].  An  important  assumption  of  these  methods 
is  that  the  incoming  signals  are  narrowband,  so  that  the  propagation  delay  of  the  signals  from  one 
subarray  to  another  can  be  approximated  by  phase  shift.  However,  if  the  FH  system  under  consid¬ 
eration  is  a  wideband  system,  then  the  inherent  frequency  variability  poses  special  difficulties  for 
signal  parameter  estimation,  due  to  the  fact  that  the  phase  shifts  among  subarrays  are  (wideband) 
frequency  dependent.  Nevertheless,  these  methods  can  still  be  applied  to  individual  hop-free  seg¬ 
ments  if  hop  timing  is  known,  since  for  a  hop-free  data  segment,  each  of  the  signals  impinging  on 
the  sensor  array  can  be  modeled  as  a  narrowband  signal. 

Here  we  extend  the  idea  of  DP-2DHR  by  coupling  the  DP  principle  with  low-rank  decompo¬ 
sition  of  data  collected  from  an  MI  sensor  array.  The  resulted  algorithm  is  termed  the  DP-TALS 
algorithm  where  TALS  stands  for  Trilinear  Alternating  Least  Squares.  Suppose  a  receiver  utilizes 
a  2-D  antenna  array,  which  is  composed  of  H  identical  subarrays  of  m  sensors  each  displaced  in 
different  directions.  Though  each  signal’s  carrier  frequencies  are  hopped  over  a  wide  frequency 
band,  between  any  two  system  wide  consecutive  hop  instants,  e.g.,  and  nl+1,  the  discrete-time 
baseband  equivalent  model  for  the  array  output  can  still  be  written  as  an  M  x  (rai+i  —  n,)  matrix 

X,  =  A,Sj  +  Wh  (10) 

where  Ai  =  [a(oti,ip i)  •  ■  ■  a(ad,  -0d)],  and  ay ,  Ty  are  azimuth  and  elevation  angles.  The  frequency- 
dependent  complex  amplitude  (due  to  path  attenuation)  is  absorbed  into  At.  Let  Jh  denote  the 
77i  x  M  selection  matrix  that  extracts  the  m  rows  corresponding  to  the  h-th  subarray,  then  it  holds 
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Figure  7:  A  rectangular  sensor  array  possessing  multiple  shift  invariances. 


that  [28] 


'  J  i  ' 

■  ' 

Yi  = 

Xi  = 

J  H 

S  +  Vt 


(11) 


where  the  m  x  d  matrix  Aj  is  the  response  of  subarray  1  (reference),  and  is  a  d  x  d  diagonal 
matrix  of  phase  shifts,  which  is  a  function  of  signal  parameters  (DOA  and  frequency)  and  the  dis¬ 
placement  of  the  h- th  subarray  relative  to  the  reference,  with  =  I.  VI'J'>  is  the  corresponding 
noise  matrix.  Define  a  H  x  cl  matrix  dy  such  that  its  h- th  row  consists  of  the  diagonal  elements  of 
then  (1 1)  can  be  rewritten  as 


Yi  =  (dy  ©  Ai)  Sj  +  Vi,  (12) 

where  0  stands  for  the  Khatri-Rao  product. 

For  given  ry  and  nl+  \ ,  the  objective  is  to  blindly  estimate  2-D  directions  or  and  0r,  as  well 
as  frequency  Up)  from  Y i,  in  (12),  for  r  —  1, ...  ,d.  A  key  observation  is  that  with  proper 
dimensioning  and  under  certain  relatively  mild  conditions,  Eqn.  (12)  is  in  fact  a  low-rank  trilinear 
(, three-way )  model  that  exhibits  strong  identifiability  properties  [23]  and  can  be  estimated  via  well- 
established  iterative  least  squares  algorithms  [24].  Low-rank  three-way  array  decomposition  is 
unique  under  a  relatively  mild  rank-like  condition  [10].  The  identifiability  of  the  model  (12)  is 
established  in  [24]. 

In  particular,  TALS  can  be  used  to  estimate  Aj,  <F , ,  and  Si  from  the  noisy  observations  Y 
The  basic  idea  of  ALS  is  to  update  matrices  one  by  one  in  an  alternating  fashion  during  each 
iteration,  conditioned  on  previously  obtained  estimates  for  the  remaining  matrices  [24].  Upon 
convergence  of  TALS,  and  Aj,  dy,  and  Si  will  be  estimated  up  to  scaling  and  common  permutation 
of  columns.  The  azimuth  and  elevation  angles  can  then  be  estimated  via  simple  division  from  Ty, 
and  the  temporal  frequencies  can  be  estimated  from  Si  via  single  1-D  harmonic  retrieval  techniques 
(e.g.,  periodogram)  or  simple  division.  Since  the  permutation  of  columns  is  common  to  all  three 
matrices,  (ar,  Vv,  u>rP^)  will  be  paired  up  automatically  by  TALS.  Notice  that  both  A;  and  Si  in 
(12)  are  Vandermonde.  This  constraint  can  be  incorporated  into  the  iteration  process  of  TALS  to 
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Figure  8:  DP-TALS:  RMSE  of  timing  and  frequency  estimation  using  a  2-D  antenna  array. 

expedite  convergence  and  improve  estimation  performance.  Joint  timing  and  frequency  estimation 
is  achieved  by  coupling  DP  with  TALS  (DP-TALS)  [13]. 

We  have  tested  the  algorithm  with  a  rectangular  array  of  size  2x6,  comprising  four  overlapping 
uniform  linear  subarrays  of  5  sensors  each.  The  spacing  A  is  half  wavelength  at  fc  =  2  GHz.  Two 
BFSK  modulated  signals  with  2-D  DOA  (elevation  ct,  azimuth  ip)  angles  of  (10°,  6°)  and  (17°, 
12°),  hop  with  different  timing.  The  hopping  bandwidth  is  80  MHz  with  1  MHz  channel  spacing. 
Both  symbol  period  and  hop  dwell  are  0.8  p,s.  The  received  signal  is  sampled  at  80  MHz  after 
down-conversion,  and  N  =  48  samples  are  collected  at  each  sensor. 

For  each  realization  of  the  Monte  Carlo  simulation,  each  of  the  two  FH-FSK  signals  hops  once 
within  the  observation  window.  Hop  timing  and  frequencies  are  randomly  generated.  The  RMSE 
versus  SNR  curves  of  the  DP-TALS  algorithm  in  Figure  8  demonstrate  that  the  algorithm  performs 
quite  well  for  a  wide  range  of  SNRs.  However,  we  note  that  the  complexity  of  DP-TALS  is  higher 
than  that  of  DP-2DHR  due  to  the  iterative  nature  of  TALS. 

B.4  An  EM  Algorithm  for  Hop  Timing  Estimation 

We  have  shown  that  the  minimizing  of  (3)  is  not  analytically  solvable,  and  one  way  to  compute 
it  is  by  dynamic  programming,  which  yield  optimal  solution,  but  comes  with  high  computational 
complexity  ( 0(KN 5))  that  may  prohibit  real-time  implementation.  The  expectation-maximization 
(EM)  principle  offers  a  second  alternative  to  solve  this  problem.  It  is  more  efficient  but  sometimes 
only  offers  sub-optimal  solutions  since  it  may  converge  to  a  local  minimum  instead  of  a  globe 
minimum.  The  EM  algorithm  was  introduced  in  the  statistics  literature  as  a  general  approach  for 
iteratively  maximizing  likelihood  functions.  It  has  found  applications  in  many  estimation  prob¬ 
lems  in  signal  processing  [19].  For  example,  the  EM  algorithm  has  been  proposed  for  sequence 
detection  and  timing  synchronization  by  several  researchers  [3,5,20,21].  [3]  discusses  sequence 
recovery  of  multiple  DS-CDMA  users  with  known  synchronization,  while  [5,20,21]  consider  joint 
sequence  detection  and  timing  estimation  for  a  single  user  transmission  of  baseband  or  direct  se¬ 
quence  CDMA  (DS-CDMA)  signals.  The  problem  studied  here  is  thus  unique  in  several  aspects: 
we  consider  joint  sequence  detection  and  timing  estimation  for  multiple  FH  spread-spectrum  trans- 
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missions;  frequency  collisions  may  be  present;  and  the  hop  pattern  is  unknown  hence  matched 
filtering  is  not  applicable. 

To  solve  (3),  let  X  represent  the  observed  but  incomplete  data,  and  (X,  n)  be  the  complete 
data.  Further  define 

0=[0T,/3Wf 

as  the  set  of  parameters  to  be  evaluated.  Using  the  EM  principle,  the  value  of  0  which  maximizes 
the  likelihood  function  p{X \<p)  corresponds  to  the  value  of  0  found  iteratively  through  steps: 


E-step:  compute 


M-step:  compute 


QU\^A  =  £[logp(X,n|0)|X,0 


ip) 


0(p+i)  _  arg  max Q  0^ 


(13) 

(14) 


the  conditional  expectation  in  the  E-step  is  with  respect  to  the  conditional  density  of  the  random 
parameter  n  given  the  X  and  assuming  that  0  =  0l!/':. 

From  Bayes  rule,  it  is  ready  to  show  that  p(X ,  n\<f> )  =  p(X\n,  0)p(n|0).  Since  the  noise 
is  white  Gaussian,  p(X\n,  0)  is  a  joint  Gaussian  pdf.  Because  n  is  independent  to  0,  p(n|0)  = 
p(n).  Thus,  E\p(n)\X,  0(p)]  does  not  depend  on  either  0  or  0(k).  After  simple  manipulation  of 
the  log  likelihood  function  and  dropping  unnecessary  terms,  we  can  rewrite  the  E-step  as: 


E-step:  compute 

where  D  =  [D0  Di 


\X-D\\2f 


X,0 


Q  (</#(p))  =  E  L 

•  •  Dk_ i],  and  Di  =  AiB.Sj. 

To  further  simplify  Q(0|0(p')),  let  us  consider: 

||X  -  D\\2F  =  tr  {(X  -  D)(X  -  D)H }  , 


ip) 


therefore 


E 


|X  —  D\\2f\X,  0^ 


=  tr  \  XX 


H 


E 


D |X,0(p)  XH 


-XE 


Dh\X,  0(p)1  +E  \dDh\X,cI)(p) 


Now  let  us  consider  one  term  of  (15).  We  have 

K- 1 


E 


D |X,0(p)  XH  =  ^E  £>j|X,  0(p) 


X 


H 


i=0 
K- 1 


Sf \x,^ 


i= 0 


X 


H 


(15) 


(16) 


since  and  do  not  depend  on  n  if  (p(p)  is  given.  Other  expectation  terms  in  (15)  can  be 
expanded  similarly  to  (16).  Based  on  (15),  as  shown  in  [13],  it  can  be  verified  that  a  simplified  EM 
algorithm  is  given  by 


16 


E-step:  for  1  <  i  <  K  —  1,  compute 


~  ip) 

n\  =  arg 


nnn 


X,  -  y 


^  tWDW/cW^r"2 

Ai-1  ~ 


(17) 


M-step:  compute 


K 


0(P+!) 

=  arg  min 

4> 


Y  x  '  A  B  s 


T\\2  I 


^  IIF 


(18) 


i=l 


where  cf/1' '  1  j  can  are  obtained  by  applying  the  MDF  algorithm  to  X;  ’s  determined  by  nil,] . 


In  summary,  given  a  received  data  block  with  multiple  hops,  the  EM  algorithm  first  takes  a  ran¬ 
dom  guess  of  ri{)}  as  the  initial  for  the  E-step,  and  compute  cj){  1 !  as  the  initial  for  the  M-step.  Then 
the  algorithm  iterates  the  following  two  steps  until  convergence:  the  expectation  step,  Eqn.  (17), 
involves  assigning  signal  segments  to  the  current  estimated  hop  frequencies  that  fits  them  best; 
The  maximization  step,  Eqn.  (18),  provides  a  new  estimate  of  hop  frequencies,  is  accomplished  by 
applying  the  MDF  algorithm  to  data  segments  according  to  the  updated  assignment.  Upon  conver¬ 
gence,  frequencies  and  complex  amplitudes  of  different  segments  pertaining  to  a  particular  signal 
can  be  associated  via  their  corresponding  DOA  parameters,  since  for  a  single  segment,  frequency 
and  DOA  parameters  pertaining  to  one  user  are  paired  up  automatically  by  the  MDF  algorithm. 
Therefore  joint  estimation  of  hop  timing  and  frequencies  of  all  users  are  obtained. 

Figure  9  depicts  the  RMSE  of  hop  timing  and  frequency  estimation  of  three  FH-FSK  signals. 
For  each  realization,  each  of  the  three  FH-FSK  signals  hops  three  times  within  the  observation 
window.  Hop  timing  is  randomly  generated,  and  frequencies  are  also  randomly  selected  from 
the  80  candidate  bins  with  the  constraint  that  there  are  always  two  collisions  randomly  located 
in  the  10  hop-free  segments  (in  average  collision  accounts  for  about  20%  time  of  the  observation 
duration).  The  RMSE  vs.  SNR  performance  shown  in  these  figures  indicates  that  the  EM  algorithm 
performs  quite  well  even  in  the  low  SNR  regime,  given  the  fact  that  the  parameters  are  estimated 
in  a  situation  where  hop  pattern,  rate,  and  timing  are  all  unknown. 


B.5  Identifiability  of  2-D  and  Multidimensional  Frequency  Estimation 

The  blind  multiuser  tracking  techniques  based  on  the  dynamic  programming  and  the  EM  principles 
described  above  require  a  coupled  2-D  frequency  estimation  step.  Therefore  the  study  of  identi¬ 
fiability  (ID)  of  2-D  frequency  estimation  is  valuable  in  determining  the  maximum  number  of 
frequencies  (thus  the  number  of  FH  signals)  that  can  be  recovered  from  a  given  data  size.  Recently 
much  progress  has  been  made  on  the  deterministic  ID  [22]  and  statistical  ID  [9, 14]  of  multidimen¬ 
sional  frequency  estimation.  For  example,  in  the  2-D  case,  the  MDF  algorithm  and  the  Unitary 
ESPRIT  algorithm  can  uniquely  estimate  up  to  0.25MiM2  frequencies  almost  surely  from  a  single 
snapshot  of  AT)  x  M2  data  mixture,  provided  that  the  2-D  frequencies  are  drawn  from  a  continuous 
distribution  (hence  the  name  of  statistical  identifiability).  In  this  project,  we  have  proved  the  most 
relaxed  statistical  identifiability  bound  to  date,  as  given  in  the  following  theorem  [12]. 


17 


Figure  9:  The  EM  algorithm:  RMSE  of  timing  and  frequency  estimation  vs.  SNR. 
Theorem  1  Given  a  sum  of  F  2-D  exponentials  as 

F 

xmi,m2  =  (19) 

/= i 

for  mi  =  1 Mi,  and  m2  =  1, . . . ,  M2,  and  without  loss  of  generality,  assume  that  > 
M2,  then  the  parameter  triples  (ujf,  Uj ,  Cf),  f  —  1,  •  •  •  ,  F,  are  almost  surely  uniquely  resolvable, 
provided  that 

F<  max  min  (2K\K2,  L  |  L  > )  (20) 

K1+L1=Mi 

K2-\-L2=M2~\-l 

where  (c Of,  Uj )  are  assumed  to  be  drawn  from  a  distribution  that  is  continuous  with  respect  to  the 
Lebesgue  measure  in  1  Er,  and  Cf  is  drawn  from  a  continuous  distribution  on  C. 

We  need  the  following  results  to  prove  Theorem  1.  Given  (19),  define  two  Vandermonde  matrices 
A  G  cMlxF  and  B  e  CM2XF,  with  generators  e]UJf  and  e]lJf ,  /  =  F ,  respectively,  then  the 

2-D  mixture  in  (19)  can  be  written  in  matrix  form  as 

X  =  AD(c)Bt,  (21) 

where  c  =  [ci,  c2, . . . ,  cp]T-  Eqn.  (21)  can  also  be  written  in  vector  form.  For  example,  let 

x  =  [xip,  Xi)2 ,  •••  ,  Xi:m2i  x2,li  •••  ,  Xm1,M2\T , 

then  it  can  be  verified  that 

x  =  (AQ  B)c,  (22) 

Proposition  1  If  we  define  a  two-dimensional  smoothing  operator  S  for  the  measurement  vector 
x  in  (22)  as 


S(x)  :=  [JiAx  •  ■  -Ji,l2x  J2, ix  ■  ■  ■  J 2,l2x  •  •  ■  JLi,l2x\, 


18 


where 


Jp 1  —  [Oa'iX(p-i)  I Ki  OjflX(Li-P)]  G  CAlxMl, 

•J©2  =  [0^2x(g-l)  /A,  0k2x(L2-?)]  £  C***, 

Jp,q  =  J?®J?\  (23) 

rmJ  1  <  p  <  Li,  1  <  q  <  L2,  and  I\,  and  L,„  for  i  —  1,2,  are  positive  integers  such  that 

K\  +  L\  =  +  1,  I\2  +  L2  =  M2  +  1.  (24) 

Then  it  can  be  verified  that 

Xs  =  Six)  =  (A(Kl)  0  B{I<2)^j  D  (c)  (A(Ll)  ©  B(L2)Y  . 

Next  we  present  the  proof  of  Theorem  1  in  the  following. 

Proof:  Given  x,  define  two  selection  matrices: 

j  1  =  0(Mi-l)xl]  ,  J 2—  [0(Mi-l)xl  -Tfcfi-l]  • 

Due  to  the  shift  invariance  property  of  Vandermonde  matrices,  we  have 

*!  =(Jx  0  Im2)x  =  ©  B)c,  (25) 

x2  ={J2  ©  IM2)x  =  (A(Ml_1)  ©  B)D[llj)c,  (26) 

ji 

where  a;  is  given  in  (22),  and  u:  :=  \ejuJI ,  e-7^2,  •  •  •  ,  eJUJF]  .  We  can  now  apply  the  2-D  smoothing 
operator  S  defined  in  Proposition  1  to  xx  and  x2.  Since  the  sizes  of  X\  and  x2  are  {M\  —  1)  x  M2, 
the  integers  in  (24)  should  be  chosen  such  that 

K\  +  L\  =  Mi,  K2  +  L2  =  M2  +  1.  (27) 

Applying  the  2-D  smoothing  operator,  we  obtain 

XhS  =  S(x1)  =  GD(c)HT , 

X  2g  =  S{x2)  =  GD{c)D{u;)Ht , 

where  G  :=  A(I<l)  ©  B[I<2\  and  H  :=  A(Ll)  ©  B(L2). 

To  further  explore  the  data  structure,  we  can  perform  the  backward  smoothing  on  the  data 
vector  x  in  (22).  Define 

y  :=  Tlx*  =  [A  0  B)c,  (28) 

where  II  is  a  backward  permutation  matrix  with  ones  on  its  antidiagonal,  and  c  =  [ci ,  c2,  ■  ■  ■  ,  cp]T, 
with 


(29) 
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Applying  the  same  technique  to  y  that  we  used  to  obtain  Xits  and  X2.s  from  x,  we  have 


Vi  =  (J i  ®  Im2)v  =  (A^1-1)  ©  B)c, 
y2  =  (J2  ®  Im2)v  =  (A(Ml_1)  ©  B)D( w)c, 
=  5(yi)  =  GD(c)HT , 

Y  2,5  =  S(y2)  =  GD(c)D(uj)Ht. 


If  we  define  the  following  matrices 


then  we  have 


Zi  :  = 


Xi,5 

.^1,5. 


Z  :  = 


'Z{ 

p 

_Z2_ 

PD(  u>) 

H 1 


(30) 


(31) 


where  the  size  of  Z  is  4K\K2  x  L1L2,  and  a  2K1K2  x  F  matrix  P  is  defined  as  P  :=  [c  c\T  0  G. 
Since  both  P  and  P  are  Khatri-Rao  products  of  Vandermonde  matrices,  invoking  the  Theorem  1 
of  [14],  when  2K\K2  >  F  and  L,  L2  >  F.  P  and  H  are  almost  surely  full  column  rank.  Hence 
Z\  and  Z2  are  of  rank  F.  The  singular  value  decomposition  of  Z  yields 


Zx 

Z2 


(32) 


where  U s  has  F  columns  which  together  span  the  column  space  of  Z.  Since  the  same  space  is 
spanned  by  the  columns  of  [(P)T  (Pfi(w))Tf ,  there  exist  an  F  x  F  nonsingular  matrix  T_1 
such  that 


U,  = 


p 

U2_ 

Pfi(w) 

-1-1 


Here  we  divide  Us  into  two  equal  size  sub-matrices  U  { .  and  U2,  each  of  size 
U\U2  satisfies: 

U\U2  =Tfi(oj)r1. 


(33) 

2KxK2  x  F.  Then 

(34) 


Assuming  that  the  elements  of  cu  are  distinct,  T  can  be  obtained  from  the  eigenvalue  decom¬ 
position  of  U\U2  up  to  column  permutation  and  scaling  ambiguity.  Suppose  that  the  eigenvalue 
decomposition  of  U\U2  gives 


T'  =TAA. 


(35) 


where  A  is  a  nonsingular  diagonal  column  scaling  matrix  and  A  is  a  permutation  matrix.  Because 
the  column  scaling  and  permutation  will  not  have  material  effect  on  the  algorithm,  for  notation 
simplicity  we  use  the  same  symbol  for  a  matrix  with  or  without  scaling  and  permutation  as  long 
as  it  is  clear  from  the  context.  Once  we  obtain  T' .  we  can  retrieve  P  and  H  up  to  common 
permutation  and  column  scaling  according  to 

P'  =  U\T'  =  PAA, 

H'  =  [p']Z^j7  =  HA  1  A. 
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(36) 

(37) 


Table  4:  Identifiability  bound  of  2-D  frequency  estimation 


M\\M2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

3 

4 

4 

6 

8 

8 

10 

12 

12 

14 

16 

4 

4 

6 

8 

9 

12 

12 

15 

16 

18 

20 

5 

6 

8 

9 

12 

12 

16 

18 

20 

21 

24 

6 

8 

9 

12 

12 

16 

18 

20 

24 

24 

28 

7 

8 

12 

12 

16 

18 

20 

24 

25 

30 

32 

8 

10 

12 

16 

18 

20 

24 

25 

30 

32 

36 

9 

12 

15 

18 

20 

24 

25 

30 

32 

36 

40 

10 

12 

16 

20 

24 

25 

30 

32 

36 

40 

42 

11 

14 

18 

21 

24 

30 

32 

36 

40 

42 

49 

12 

16 

20 

24 

28 

32 

36 

40 

42 

49 

50 

In  (37)  we  have  used  the  fact  that  A-1  =  A1 .  The  permutation  is  not  an  issue  here,  however, 
because  c Of  and  Uf  appear  in  the  same  column  of  P',  as  well  as  in  H\  thus  are  automatically  paired; 
and  arbitrary  nonzero  column  scaling  is  immaterial,  because  the  frequencies  can  be  obtained  by 
dividing  suitably  chosen  elements  of  the  aforementioned  column.  For  example,  suppose  eJUJf  and 
e]1Jf  appear  in  the  /- th  column  of  P',  then  e]LOf  can  be  retrieved  by  any  one  of  the  following 
equations  in  the  noiseless  case: 

p' 

ej“f  =  ,  n  =  K2{h  -  1)  +  k2,  (38) 

Pn-K2J 

where  k\  —  2,  •  •  •  ,  Ku  k2  =  1,  •  •  •  ,  K2,  and  similarly,  ev'!  can  be  retrieved  by  any  one  of  the 
following: 

ejuf=Pn!^  n  =  K2(k1~l)  +  k2,  (39) 

Pn-l,f 

where  k\  =  1,  ■  ■  ■  ,  K\ .  and  k2  —  2,  •  •  •  ,  K> .  Notice  that  and  e3,/{  are  automatically  paired 
since  they  are  obtained  from  the  same  column  of  P' . 

Hence  we  have  shown  that  the  2-D  frequencies  can  be  uniquely  recovered  almost  surely,  pro¬ 
vided  that  F  <  2K\K2  and  F  <  L2L2,  where  Kt  and  L,  are  positive  integers  subject  to 

K\  +  L\  =  .Mi,  K2  +  L2  =  M2  +  1. 


Therefore  Theorem  1  is  proved.  ■ 

It  is  difficult  to  obtain  the  exact  solution  of  the  integer  optimization  problem  in  (20).  Here  we 
have  listed  the  maximum  number  of  2-D  frequencies  that  can  be  uniquely  identified  for  certain 
data  sizes  in  Table  4.  We  have  also  found  the  lower  and  upper  bounds  of  the  maximum  number 
identifiable  2-D  frequencies  as  given  in  the  following  proposition  [12]. 
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Proposition  2  The  maximum  number  of  2 -D  frequencies  given  in  Theorem  1  is  bounded  by 


nun 


in  {2  [(V2  -  1  )Mi]  [(\/2  -  1  )(M2  +  1)]  ,  [(2  -  [(2  -  y/2 )(M2  +  1)]  } 

min  (2K1K2,  LiLf)  <  |_0.343Mi(M2  +  1)J  . 


<  max 

A'i+Li=Mi 
i^2+^2—  -M2  +  I 


Theorem  1  shows  significant  improvement  on  the  identifiability  bound  of  2-D  frequency  es¬ 
timation  over  existing  algebraic  algorithms.  Previously  the  most  relaxed  statistical  ID  bound  is 
achieved  by  the  MDF  algorithm  as  shown  in  [15].  Using  Theorem  1  of  [14],  we  are  also  able  to 
obtain  the  statistical  ID  bounds  of  the  Unitary  ESPRIT  algorithm  [6]  and  the  MEMP  algorithm  [8]. 
Figure  10  plots  a  comparison  of  the  statistical  ID  bounds  of  different  algebraic  algorithms  for  2- 
D  frequency  estimation  in  the  absence  of  noise.  The  Unitary  ESPRIT  algorithm  and  the  MDF 
algorithm  have  the  same  ID  bound,  which  is 

F  <  \M1/ 2]  \M2/ 2]  . 


The  statistical  ID  bound  of  the  MEMP  algorithm  is  slightly  smaller  than  those  of  the  MDF  algo¬ 
rithm  and  the  Unitary  ESPRIT  algorithm  because  no  back  ward-  forward  smoothing  is  used  in  the 
MEMP  algorithm.  Note  that  the  deterministic  ID  bound  of  the  MEMP  algorithm  is 

F  <  min  (Mi/2,  M2/2) 


as  shown  in  [8]  and  [33].  It  can  be  seen  from  Fig.  10  that  Theorem  1  offers  a  significantly  improved 
ID  bound  over  existing  results.  We  note  that  an  algebraic  algorithm  for  2-D  frequency  estimation 
can  be  obtained  from  the  constructive  proof  of  Theorem  1,  which  may  be  used  to  replace  the 
aforementioned  MDF  algorithm  [12]. 

The  2-D  identifiability  results  can  be  generalized  to  the  case  of  iV-D  frequency  estimation. 


Theorem  2  Given  a  sum  of  F  N-D  exponentials 

F  N 

xm, . „„  -  y  .  n  (40) 

/= 1  n=l 

formn  —  1, . . . ,  Mn,  n  —  1, . . . ,  N,  and  without  loss  of  generality,  assume  that  Mi  =  max{Mn,  n  = 

(N  N  \ 

2  IT  Kn,  IT  Ln  ) ,  (41) 

11  11 

n  1  n  1  / 

2<n<N 

and  the  distributions  used  to  draw  the  N F  frequencies  and  F  amplitudes  are  continuous  with 
respect  to  Lebesgue  measure  in  I  l  w  and  C,  respectively,  then  the  parameter  (N  +  1) -tuples 
(c Ufti, . . .  Cf),  f  —  1, . . . ,  F,  are  almost  surely  uniquely  resolvable. 
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Figure  10:  Comparison  of  identifiability  results  for  2-D  frequency  estimation  {Mx  =  M2). 

B.6  Conclusion 

In  this  project  we  have  developed  a  signal  processing  scheme  for  blind  tracking  of  multiple  FH 
signals  over  multipath  channels.  This  technique  is  based  on  the  principle  of  dynamic  programming 
and  expectation-maximization,  coupled  with  multidimensional  harmonic  and  low-rank  analysis. 
Numerical  simulations  demonstrate  its  capability  of  joint  estimation  of  hop  timing,  frequency,  and 
DOA  of  multiple  FH  signals  in  the  presence  of  frequency  collisions,  without  the  knowledge  of 
signal  hop  patterns.  The  significance  of  this  project  in  basic  research  lies  in  innovative  methods  for 
signal  detection  and  jammer  localization  in  frequency  hopping  communications,  and  fundamental 
understanding  of  the  identifiability  of  multidimensional  frequency  estimation.  Furthermore,  it 
has  yielded  practical  algorithms  that  are  applicable  in  the  presence  of  modulation  uncertainly  and 
unknown  channels. 

B.7  Bibliography 

[1]  L.  Aydin  and  A.  Polydoros,  “Joint  hop-timing  estimation  for  FH  signals  using  a  coarsely 
channelized  receiver,”  Proc.  Military  Communications  Conference,  vol.  2,  pp.  769-773,  San 
Diego,  CA,  Nov.  1995. 

[2]  L.  Aydin  and  A.  Polydoros,  “Hop-timing  estimation  for  FH  signals  using  a  coarsely  channel¬ 
ized  receiver,”  IEEE  Trans.  Communications,  44(4):5 16—526,  Apr.  1996. 

[3]  M.  J.  Borran  and  M.  Nasiri-Kenari,  “An  efficient  detection  technique  for  synchronous  CDMA 
communication  systems  based  on  the  expectation  maximization  algorithm,”  IEEE  Trans. 
Vehicular  Technology ,  Vol.  49,  No.  5,  pp.  1663-1668,  Sept.  2000. 

[4]  U.-C.G.  Fiebig,  “An  algorithm  for  joint  detection  in  fast  frequency  hopping  system,”  Proc. 
ICC  1996 ,  vol.  1,  pp.  90-95. 

[5]  C.  N.  Georghiades  and  D.  L.  Snyder,  “The  expectation-maximization  algorithm  for  symbol 
unsynchronized  sequence  detection,”  IEEE  Trans.  Communications,  Vol.  39,  No.  1,  pp.  54- 
61,  Jan.  1991. 


23 


[6]  M.  Haardt  and  J.  A.  Nossek,  “Simultaneous  Schur  decomposition  of  several  nonsymmet- 
ric  matrices  to  achieve  automatic  pairing  in  multidimensional  harmonic  retrieval  problems,” 
IEEE  Trans.  Signal  Processing,  vol.  46,  no.  1,  pp.  161-169,  Jan.  1998. 

[7]  P.  Hande,  L.  Tong,  and  A.  Swami,  “Multipath  delay  estimation  for  frequency  hopping  sys¬ 
tems,”  Journal  of  VLSI  Signal  Processing,  vol.  30,  no.  1-3,  pp.  163-178,  2002. 

[8]  Y.  Hua,  “Estimating  two-dimensional  frequencies  by  matrix  enhancement  and  matrix  pencil,” 
IEEE  Trans.  Signal  Processing,  40(9):2267-2280,  Sep.  1992. 

[9]  T.  Jiang,  N.  D.  Sidiropoulos,  and  J.  M.  F.  ten  Berge,  “Almost-sure  identifiability  of  mul¬ 
tidimensional  harmonic  retrieval,”  IEEE  Trans.  Signal  Processing,  49(9):  1849-1859,  Sep. 
2001 

[10]  J.  B.  Kruskal,  “Three-way  arrays:  rank  and  uniqueness  of  trilinear  decompositions,  with 
application  to  arithmetic  complexity  and  statistics,”  Linear  Algebra  and  Its  Applications, 
18:95-138,  1977. 

[11]  J.  Li,  P.  Stoica,  and  D.  Zheng,  “An  efficient  algorithm  for  two-dimensional  frequency  esti¬ 
mation,”  Multidimensional  Systems  and  Signal  Processing,  7(2):  15 1-178,  Apr.  1996. 

[12]  J.  Liu  and  X.  Liu,  “An  eigenvector-based  algebraic  approach  for  two-dimensional  frequency 
estimation  with  improved  identifiability,”  IEEE  6th  SPAWC,  pp.  655-659,  New  York  City, 
NY,  Jun.  2005. 

[13]  X.  Liu,  Z.  Li,  and  X.  Ma,  “An  EM  algorithm  for  hop  timing  estimation  in  FH  networks  with 
frequency  collisions,”  Proc.  Military  Communications  Conference,  pp.  1635-1639,  Mon¬ 
terey,  CA,  Oct.  2004. 

[14]  X.  Liu  and  N.  D.  Sidiropoulos,  “Almost  sure  identifiability  of  mulitidimensional  constant 
modulus  harmonic  retrieval,”  IEEE  Trans.  Signed  Processing,  Vol.  50,  no.  9,  pp.  2366-2368, 
Sep.  2002. 

[15]  X.  Liu  and  N.  D.  Sidiropoulos,  “On  constant  modulus  multidimensional  harmonic  retrieval,” 
Proc.  ICASSP,  Orlando,  FL,  May  2002,  pp.  2977-2980. 

[16]  X.  Liu,  N.  D.  Sidiropoulos,  and  A.  Swami,  “Joint  hop  timing  and  frequency  estimation  for 
collision  resolution  in  frequency  hopped  networks,”  IEEE  Trans.  Wireless  Communications, 
vol.  4,  no.  6,  to  appear,  Nov.  2005. 

[17]  X.  Liu,  N.  D.  Sidiropoulos,  and  A.  Swami,  “Blind  high  resolution  localization  and  tracking 
of  multiple  frequency  hopped  signals,”  IEEE  Trans.  Signed  Processing,  vol.  50,  no.  4,  pp. 
889-901,  Apr.  2002. 

[18]  T.  Mabuchi,  R.  Kohno,  and  H.  Imai,  “Multiuser  Detection  Scheme  Based  on  Canceling 
Cochannel  Interference  for  MFSK/FH-SSMA  System,”  IEEE  JSAC,  12(4):593-604,  May 
1994. 

[19]  T.  K.  Moon,  “The  expectation-maximization  algorithm,”  IEEE  Signed  Processing  Magazine, 
Vol.  13,  No.  6,  pp.  47-60,  Nov.  1996. 


24 


[20]  C.  R.  Nassar  and  M.  R.  Soleymani,  “Joint  sequence  detection  and  phase  estimation  using  the 
EM  algorithm,”  Proc.  Canadian  Conference  on  Electrical  and  Computer  Engineering ,  vol. 
1,  pp.  296-299,  Sept.  1994. 

[21]  I.  Sharfer  and  A.  O.  Hero,  “Spread  spectrum  sequence  estimation  and  bit  synchronization 
using  an  EM-type  algorithm,”  Proc.  ICASSP ,  vol.  3,  pp.  1864-1867,  May  1995. 

[22]  N.  D.  Sidiropoulos,  “Generalizing  Caratheodory’s  uniqueness  of  harmonic  parameterization 
to  N  dimensions,”  IEEE  Trans.  Information  Theory,  vol.  47,  no.  4,  pp.  1687-1690,  May 
2001. 

[23]  N.  D.  Sidiropoulos  and  R.  Bro,  “On  the  uniqueness  of  multilinear  decomposition  of  N-way 
arrays,”  J.  of  Chemome tries,  14(3):229— 239,  May  2000. 

[24]  N.  D.  Sidiropoulos,  R.  Bro,  and  G.  B.  Giannakis,  “Parallel  factor  analysis  in  sensor  array 
processing,”  IEEE  Trans.  Signed  Processing,  48(8):2377— 2388,  Aug.  2000. 

[25]  M.  K.  Simon,  U.  Cheng,  L.  Aydin,  A.  Polydoros,  and  B.  K.  Levitt,  “Hop  timing  estimation  for 
noncoherent  frequency-hopped  M-FSK  intercept  receivers,”  IEEE  Trans.  Communications, 
vol.  43,  no.  2/3/4,  pp.  1144-1154,  Feb./Mar./Apr.  1995. 

[26]  A.  Souloumiac,  “Blind  Source  Detection  and  Separation  Using  Second-Order  Non- 
stationarity,”  Proc.  ICASSP,  vol.  3,  pp.  1912-1915,  1995. 

[27]  A.  L.  Swindlehurst,  B.  Ottersten,  R.  Roy,  and  T.  Kailath,  “Multiple  invariance  ESPRIT,” 
IEEE  Trans.  Signed  Processing,  vol.  40,  no.  4,  pp.  867-881,  Apr.  1992. 

[28]  A.  L.  Swindlehurst  and  T.  Kailath,  “Azimuth/elevation  direction  finding  using  regular  array 
geometries,”  IEEE  Trans.  Aerospace  and  Electronic  Systems,  vol.  29,  No.  1,  pp.  145-156, 
Jan.  1993. 

[29]  D.  Torrieri,  “An  anticipative  adaptive  array  for  frequency-hopping  communications,”  IEEE 
Trans.  AES,  24:449-456,  Jul.  1988. 

[30]  D.  Torrieri  and  K.  Bakhru,  “Frequency  compensation  in  an  adaptive  antenna  system  for 
frequency-hopping  communications,”  IEEE  Trans.  AES,  23:448-467,  Jul.  1987. 

[31]  D.  Torrieri  and  K.  Bakhru,  “Performance  of  recursive  suppression  algorithm  in  nonstationary 
environments,”  IEEE  Trans.  AES,  43:299-307,  Mar.  1995. 

[32]  K.  T.  Wong,  “Blind  beamforming/geolocation  for  wideband-FFHs  with  unknown  hop- 
sequences,”  IEEE  Trans.  AES,  37(1):65— 76,  Jan.  2001. 

[33]  H.  Yang  and  Y.  Hua,  “On  rank  of  block  Hankel  matrix  for  2-D  frequency  detection  and 
estimation,”  IEEE  Trans.  Signal  Processing,  vol.  44,  no.  4,  pp.  1046-1048,  Apr.  1996. 

[34]  M.  D.  Zoltowski  and  K.  T.  Wong,  “ESPRIT-based  2-D  direction  finding  with  a  sparse  uniform 
array  of  electromagnetic  vector  sensors,”  IEEE  Trans.  Signal  Processing,  vol.  48,  no.  8,  pp. 
2195-2204,  Aug.  2000. 


25 


C  Publications  from  this  Project 

(a)  Papers  published  in  peer-reviewed  journals: 

1.  X.  Liu,  N.  D.  Sidiropoulos,  and  A.  Swami,  “Joint  hop  timing  and  frequency  estimation 
for  collision  resolution  in  frequency  hopped  networks,”  IEEE  Trans.  Wireless  Communi¬ 
cations,  vol.  4,  no.  6,  to  appear,  Nov.  2005  (submitted  Apr.  2004). 

(b)  Book  chapters  published: 

1.  X.  Liu,  N.  D.  Sidiropoulos,  and  T.  Jiang,  “Multidimensional  harmonic  retrieval  with  ap¬ 
plications  in  MIMO  wireless  channel  sounding,”  in  A.  B.  Gershman  and  N.  D.  Sidiropou¬ 
los,  editors,  Space-Time  Processing  for  MIMO  Communications,  Wiley,  pp.  41-75,  2005. 

2.  X.  Liu  and  N.  D.  Sidiropoulos,  “PARAFAC  techniques  for  high-resolution  array  process¬ 
ing,”  in  Y.  Hua,  A.  Gershman,  and  Q.  Cheng,  editors,  High  Resolution  and  Robust  Signal 
Processing,  Marcel  Dekker  Inc.,  New  York,  pp.  111-150,  2003. 

(c)  Papers  published  in  conference  proceedings: 

1.  J.  Liu  and  X.  Liu,  “An  eigenvector-based  algebraic  approach  for  two-dimensional  fre¬ 
quency  estimation  with  improved  identihability,”  IEEE  6th  Workshop  on  Signal  Process¬ 
ing  Advances  in  Wireless  Communications,  pp.  655-659,  New  York  City,  NY,  Jun.  2005. 

2.  J.  Li,  X.  Liu,  and  A.  Swami,  “Detection  of  model  order  variation  in  frequency  hopping 
system  with  bandwidth  mismatch,”  IEEE  6th  Workshop  on  Signal  Processing  Advances 
in  Wireless  Communications,  pp.  680-684,  New  York  City,  NY,  Jun.  2005. 

3.  X.  Ma,  S.  Zhou  and  X.  Liu,  “A  novel  view  on  modulus-preserving  rate-one  space  time 
block  codes,”  Proc.  of  the  2005  International  Conf  on  Wireless  Networks,  Communica¬ 
tions,  and  Mobile  Computing,  Maui,  HI,  Jun.  2005. 

4.  X.  Liu,  Z.  Li,  and  X.  Ma,  “An  EM  algorithm  for  hop  timing  estimation  in  FH  networks 
with  frequency  collisions,”  Proc.  Military  Communications  Conference,  pp.  1635-1639, 
Monterey,  CA,  Oct.  2004. 

5.  X.  Liu,  Z.  Li,  N.  D.  Sidiropoulos,  and  A.  Swami,  “Joint  signal  parameter  estimation 
of  wideband  frequency  hopped  transmissions  using  2-D  antenna  arrays,”  Proc.  Signed 
Processing  Advances  in  Wireless  Communications,  pp.  624-628,  Rome,  Italy,  Jun.  2003. 

(d)  Manuscripts  submitted,  but  not  published: 

1.  J.  Liu  and  X.  Liu,  “An  eigenvector  based  algebraic  approach  for  multidimensional  fre¬ 
quency  estimation  with  improved  identihability,”  IEEE  Trans.  Signed  Processing,  sub¬ 
mitted  May  2005,  revised  Oct.  2005. 

2.  X.  Liu,  J.  Li,  and  X.  Ma,  “An  EM  algorithm  for  blind  hop  timing  estimation  of  multiple 
frequency  hopped  signals  with  bandwidth  mismatch,”  IEEE  Trans.  Vehicular  Technol¬ 
ogy,  submitted  Jun.  2005. 


26 


(e)  Papers  presented  at  meetings,  but  not  published  in  conference  proceedings: 

1.  X.  Liu,  “Code-blind  reception  of  frequency  hopped  signals  over  multipath  fading  chan¬ 
nels,”  at  the  International  Conference  on  Acoustics,  Speech,  and  Signal  Processing ,  Mon¬ 
treal,  Quebec,  Canada,  May  2004. 

2.  X.  Liu,  “On  constant  modulus  multidimensional  harmonic  retrieval,”  at  the  Trilinear 
Methods  in  Chemistry  and  Psychology  Conference,  Lexington,  KY,  Jun.  2003. 

(f)  Technical  reports  submitted  to  ARO: 

1.  X.  Liu,  “Signal  detection  and  jammer  localization  in  multipath  channels  for  frequency 
hopping  communications,”  Technical  report  submitted  to  ARO,  Oct.  2005. 


D  Personnel  Supported 

1.  Jingli  Li  (graduate  student) 

2.  Jun  Liu  (graduate  student) 

3.  Xiangqian  Liu  (principle  investigator) 

E  Inventions  and  Patents 

None. 


27 


